1. Field of the Invention
The present invention relates to a method of measuring strain in a target body, using the transmission, reception, processing and normalization of ultrasound signals.
2. Description of the Prior Art
Imaging of elastic parameters of soft tissue has developed into a new tool for diagnosis of disease. Current estimators of tissue motion, include a time-domain cross-correlation based speckle tracking algorithm, and a Fourier based speckle phase-tracking technique. These techniques are coherent estimation techniques, i.e., these methods are sensitive to phase variations. The coherent estimation techniques generally have the advantage of being highly precise. Strain Filter (SF) analysis has shown, however, that they are not very robust in the presence of even a small amount of decorrelation between the pre- and post-compression signals. The term xe2x80x98robustnessxe2x80x99 has been used in statistical analysis to denote the good performance of statistical tests, i.e., the homogeneity of the variance calculation, even if the data deviates from theoretical requirements. By equivalence in elastography, xe2x80x3robustnessxe2x80x2 denotes the consistently good performance of the estimator even at high decorrelation noise, i.e., keeping the variance of estimation at a relatively constant and low level at a large range of noise levels.
The term xe2x80x9cdecorrelationxe2x80x9d as used herein, is defined as the loss of full correlation between the pre- and postcompressed windowed signal segments. Therefore, decorrelation may be encountered due to many sources, such as intrawindow axial motion undesired lateral or elevational motion, jitter (i.e., any cause of misregistration between the pre- and post-compressed A-line segments), unstable mechanical setup, etc. The main idea in this study is to introduce a new estimator that is more immune to decorrelation compared to other estimators.
The tissue strain estimator is a spectral estimator that estimates strain directly. Since the proposed estimator uses the power spectrum, it is incoherent, i.e., it does not use the phase of the signal. Previously reported incoherent methods include optical flow speckle tracking, envelope cross-correlation, and spectral chirp z-transform techniques. Generally, incoherent methods may be less precise but significantly more robust. For example, we have demonstrated this property for the case of time-delay estimation using the envelope of echo-signals. This may be a significant advantage where elastography is to be practiced in situations involving (1) undesired scanning motion, such as the case of using an unstable handheld transducer and/or (2) undesired tissue motion, such as abdominal or intravascular elastography. This property of the estimator is demonstrated later in this paper in the simulation results section through testing of its immunity to noise caused by jitter.
The main idea behind a spectral strain approach is based on the Fourier scaling property, which implies that a compression or expansion of the time-domain signal should lead to an expansion or compression of its power spectrum, respectively. One of the most well known and thoroughly studied spectral motion effects is the Doppler shift, which typically links the frequency shift to the scatterer velocity between emissions. Velocities towards the transducer result in a positive frequency shift, while the opposite is true for scatterers that move away from the transducer. However, since the scatterers within a given resolution length do not move at the same velocities, a spectrum of Doppler frequencies is observed. Therefore, initially in ultrasound, the methods of velocity estimation for the measurement of blood flow mainly operated in the frequency domain, otherwise known as spectrum analysis techniques and measured the mean velocity of scatterers across the vessel lumen (indicative of the volumetric flow rate) by estimating the mean frequency of the power spectrum. Despite the success of these techniques even in in vivo vessels, detection of the Doppler frequency shift, which is typically on the order of 1 kHz, is not possible for pulsed instruments, since the downshift in frequency due to attenuation (on the order of 10-100 kHz) is expected to dominate over the Doppler shift. Since in elastography the pre- and postcompressed segments are approximately identical depths, the attenuation effect on the two spectra is assumed to be identical and cancelled out when the two spectra are compared.
Strain estimation using spectral methods depends on the subsequent change in the scatterer statistics. Spectral methods typically link one or more signal parameters to the change in mean scatterer spacing. One prior art relates the relative change in the mean scatterer spacing to the strain incurred during a cardiac cycle. This method assumes the presence of underlying scatterer periodicities. Despite the fact that this has also been demonstrated to work in in vivo intravascular applications, the main assumptions of regular spacing or periodicities may not hold for most tissues. In contrast, as shown in the theory section, the spectral methods mentioned in this paper make no assumptions regarding the composition of the tissue scatterers.
Typically in elastography time-domain techniques are used that involve the computation of the time-delay to estimate the displacement following an applied compression, and the estimation of strain by applying gradient operations on the previously obtained time-delay estimates. As mentioned earlier, an important advantage associated with these spectral methods as well as other estimators, such as the adaptive stretching estimator, is that they can be used to estimate strain directly; i.e., without involving the use of noisy gradient operators. In the latter case, the gradient operation introduces additional amplification of the noise into the strain estimation process, thus degrading the strain estimates. Furthermore, similar to the adaptive stretching estimator, only one estimation window is needed, for both the magnitude and the sign of the strain to be estimated.
As shown later in the theory section, spectral estimators can be divided into two main groups: a) the spectral shift methods and b) the spectral bandwidth methods. Despite the fact that we develop expressions that show direct strain estimation in both cases, in this paper we focus primarily on a spectral shift method; we estimate the relative shift in the spectral centroid caused by compressive or tensile tissue strain. Therefore, throughout this paper this new estimator is referred to as the xe2x80x9ccentroid strain estimatorxe2x80x9d, xe2x80x9ccentroid estimatorxe2x80x9d or xe2x80x9ccentroid methodxe2x80x9d. Current investigations that deal with the development of alternative shift estimators as well as bandwidth estimators will not be reported in this study.
The spectral centroid has been widely used in estimating the Doppler shift, attenuation and backscattering. The theory underlying the use of centroid strain estimators is presented in the next section. One-dimensional (1-D), instead of two-dimensional (2-D), motion simulations are used in order to more accurately study the performance of the estimator, i.e., independent of the effect of signal decorrelation in two dimensions that complicates the measurements. Simulation results in 1-D illustrate the insensitivity of the centroid strain estimator to signal decorrelation effects. It is important to note that, as mentioned earlier, decorrelation can be due to several sources. For the purpose of this paper, we consider solely the axial decorrelation effect in this 1-D model. We, thus, assume that the robustness demonstrated by the spectral estimator vis-à-vis this effect is a more general property that can be further applied at other decorrelation scenarios. For example, it is shown in the results section how the spectral method is indeed more immune to jitter, another source of decorrelation. The elastograms obtained using these simulations as well as phantom experiments illustrate the robustness of the spectral centroid strain estimator. The properties of the new estimator are discussed and summarized in the Conclusion section.
In this section, we show analytically that for Gaussian echo spectra, the relative spectral shift is a direct measure of tissue strain. We also show the relative bandwidth variation can also be used as a direct strain estimator.
Signal and noise model.
The pre- and post-compression echo signals are given as follows:
r1(z)=h(z)*e(z)+n1(z)xe2x80x83xe2x80x83(1) 
r2(z)=h(z)*e(az)+n2(z)xe2x80x83xe2x80x83(2)
where z is a spatial variable, r1(z) and r1(z) are the received RF signals before and after compression, respectively, h(z) is the impulse response of the ultrasound system or point-spread function (PSF), e(z) is the scattering function, n1(z) and n2(z) are independent zero-mean white noise sources and xcex1 is the compression coefficient (or, strain factor) linked to strain s through                     a        =                              1                          1              -              s                                ≅                      1            +            s                                              (        3        )            
The approximation holds for s less than  less than 1, where the strain s for a one-dimensional homogeneous target is typically defined in mechanics by                     s        =                                            L              0                        -            L                                L            0                                              (        4        )            
where L0 and L are the pre- and postcompressed axial dimensions of the target. From Eq. (4) the reader should note that positive strain denotes compression (and a greater than 1) while negative strain denotes tension (and a less than 1). The reader should note that throughout this paper the subscripts 1 and 2 denote pre- and postcompression parameters, respectively.
Assuming that h(z) and e(z) in Eqs. (1) and (2) can be described by their autocorrelation functions that may be modeled by modulated Gaussian functions, we obtain                                           h            ⁡                          (              z              )                                =                                    1                                                                    2                    ⁢                    π                                                  ⁢                                  L                  h                                                      ⁢                          exp              ⁡                              (                                                                            -                                              z                        2                                                              /                    2                                    ⁢                                      L                    h                    2                                                  )                                      ⁢            sin            ⁢                          xe2x80x83                        ⁢                          (                                                k                  h                                ⁢                z                            )                                      ⁢                  
                ⁢        and                            (        5        )                                          e          ⁡                      (            z            )                          =                              1                                                            2                  ⁢                  π                                            ⁢                              L                e                                              ⁢                      exp            ⁡                          (                                                                    -                                          z                      2                                                        /                  2                                ⁢                                  L                  e                  2                                            )                                ⁢          cos          ⁢                      xe2x80x83                    ⁢                      (                                          k                e                            ⁢              z                        )                                              (        6        )            
where Lh and Le are the resolution lengths of the PSF and of the scattering function, respectively, kh is the central spatial frequency of the PSF, and ke is the central spatial frequency of the scattering function.
The one-sided power spectra of the pre- and post-compression RF signals (positive frequencies) are given respectively by                                           R            1                    ⁡                      (            k            )                          =                                            1              4                        ⁢                          exp              ⁡                              (                                  -                                                            1                      2                                        ⁡                                          [                                                                                                                                  (                                                              k                                -                                                                  k                                  h                                                                                            )                                                        2                                                    ⁢                                                      L                            h                            2                                                                          +                                                                                                            (                                                              k                                -                                                                  k                                  e                                                                                            )                                                        2                                                    ⁢                                                      L                            e                            2                                                                                              ]                                                                      )                                              +                                    N              l                        ⁡                          (              k              )                                                          (        7        )                                                      R            2                    ⁡                      (            k            )                          =                                            1                              4                ⁢                a                                      ⁢                          exp              ⁡                              (                                  -                                                            1                      2                                        ⁡                                          [                                                                                                                                  (                                                              k                                -                                                                  k                                  h                                                                                            )                                                        2                                                    ⁢                                                      L                            h                            2                                                                          +                                                                                                            (                                                              k                                -                                                                  ak                                  e                                                                                            )                                                        2                                                    ⁢                                                                                    L                              e                              2                                                                                      a                              2                                                                                                                          ]                                                                      )                                              +                                    N              2                        ⁡                          (              k              )                                                          (        8        )            
where N1(k) and N2(k) are independent power spectra of zero-mean white noise processes, i.e.,                     ⟨        N            1              (        k        )              ⟩    =                              ⟨          N                2                  (          k          )                    ⟩        =    0.  
A brief observation of Eqs. (7) and (8) reveals the centroid shift in the scattering spectrum resulting from the compression. In other words, if fe1 and fe2 are the center frequencies of the scattering spectrum before and after compression, respectively, and assuming that the speed of sound in the tissue c remains constant, from Eqs. (7) and (8) we have                                                                                           f                  e2                                -                                  f                  e1                                            =                                                c                                      2                    ⁢                    π                                                  ⁢                                  (                                                            ak                      e                                        -                                          k                      e                                                        )                                                                                                                        =                                                      (                                          a                      -                      1                                        )                                    ⁢                                      f                    e1                                                              ,                                                          (        9        )            
or, from Eq. (3)                                                         f              e2                        -                          f              e1                                            f            e1                          ≅        s                            (        10        )            
So, the relative centroid shift in the scattering function spectrum constitutes a direct strain estimator. A similar result is found later (Eq. 16) using the spectrum of the received signal. In the section below, we use Eq. (10) as a guide in the formulation of the new estimator.
Effect of strain on the spectrum of the received signal
The centroid of the power spectrum of the received signal is defined as follows                               f          c                =                              c                          2              ⁢              π                                ⁢                                                                      ∫                                      -                    ∞                                    ∞                                ⁢                                                                            kR                      1                                        ⁡                                          (                      k                      )                                                        ⁢                                      ⅆ                    k                                                                                                ∫                                      -                    ∞                                    ∞                                ⁢                                                                            R                      1                                        ⁡                                          (                      k                      )                                                        ⁢                                      ⅆ                    k                                                                        .                                              (        11        )            
The centroid estimate for the precompression power spectrum is given by:                               f          e1                =                              c                          2              ⁢              π                                ⁢                                                                      k                  h                                ⁢                                  L                  h                  2                                            +                                                k                  e                                ⁢                                  L                  e                  2                                                                                    L                h                2                            +                              L                e                2                                                                        (        12        )            
In a similar manner we can derive the expression for the centroid of the post-compression power spectrum by replacing ke, and Le, in Eq. (9) by their corresponding parameters in the post-compression power spectrum, i.e., ake and       L    e    a
(as indicated from Eq. 7) to obtain:                               f          c2                =                              c                          2              ⁢              π                                ⁢                                                                      k                  h                                ⁢                                  L                  h                  2                                            +                                                k                  e                                ⁢                                                      L                    e                    2                                    a                                                                                    L                h                2                            +                                                L                  e                  2                                                  a                  2                                                                                        (        13        )            
Note that the PSF parameters remain unchanged. Since both centroids depend on the center frequencies and bandwidths of the scattering function and the PSF and by consulting Eq. (10), we normalize this effect by using the following ratio as a candidate strain estimator:                                                         f              c2                        ⁢                          xe2x80x83                        -                          xe2x80x83                        ⁢                          f              c1                                            f            c1                          ⁢                  xe2x80x83                =                  xe2x80x83                ⁢                                                                                                  k                    h                                    ⁢                                      xe2x80x83                                    ⁢                                      a                    2                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    h                    2                                                  ⁢                                  xe2x80x83                                +                                  xe2x80x83                                ⁢                                                      ak                    e                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    e                    2                                                                                                                    a                    2                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    h                    2                                                  ⁢                                  xe2x80x83                                +                                  xe2x80x83                                ⁢                                  L                  e                  2                                                      ⁢                          xe2x80x83                        -                          xe2x80x83                        ⁢                                                                                k                    h                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    h                    2                                                  ⁢                                  xe2x80x83                                +                                  xe2x80x83                                ⁢                                                      k                    e                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    e                    2                                                                                                L                  h                  2                                ⁢                                  xe2x80x83                                +                                  xe2x80x83                                ⁢                                  L                  e                  2                                                                                                                          k                  h                                ⁢                                  xe2x80x83                                ⁢                                  L                  h                  2                                            ⁢                              xe2x80x83                            +                              xe2x80x83                            ⁢                                                k                  e                                ⁢                                  xe2x80x83                                ⁢                                  L                  e                  2                                                                                    L                h                2                            ⁢                              xe2x80x83                            +                              xe2x80x83                            ⁢                              L                e                2                                                                        (        14        )            
The parameters Le and Lh are related to Be and Bh, the equivalent noise spectral bandwidths for the scattering and PSF spectra, through             B      h        =                            1                      2            ⁢                          π                        ⁢                          L              h                                      ⁢                  xe2x80x83                ⁢        and        ⁢                  xe2x80x83                ⁢                  B          e                    =              1                  2          ⁢                      π                    ⁢                      L            e                                ,
respectively.
However, the PSF bandwidth is typically much smaller than the bandwidth of the scattering function, i.e., Be greater than  greater than Bh; thereby, Le less than  less than Lh and, therefore, a2L2h greater than  greater than L2h. After cancellation of common terms in the numerator and denominator of Eq. (14), we obtain                                                                                                               f                    c2                                    -                                      f                    c1                                                                    f                  c1                                            ≅                              xe2x80x83                            ⁢                                                                    (                                          a                      -                      I                                        )                                    ⁢                                                                                    k                        e                                            ⁢                                              L                        e                        2                                                                                                            L                        h                        2                                            +                                              L                        e                        2                                                                                                                                                                                k                        h                                            ⁢                                              L                        h                        2                                                              +                                                                  k                        e                                            ⁢                                              L                        e                        2                                                                                                                        L                      h                      2                                        +                                          L                      e                      2                                                                                                                                              =                              xe2x80x83                            ⁢                                                (                                      a                    -                    1                                    )                                                                                                                    k                        h                                            ⁢                                              L                        h                        2                                                                                                            k                        e                                            ⁢                                              L                        e                        2                                                                              +                  1                                                                                        (        15        )            
or, from the small strain approximation case of Eq. (3) (i.e., in mathematical terms, for strains less than 10%),                                                         f              c2                        -                          f              c1                                            f            c1                          ≅        As                            (        16        )            
where A is given by                               A          ⁢                      xe2x80x83                    =                      xe2x80x83                    ⁢                      1                                                                                k                    h                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    h                    2                                                                                        k                    e                                    ⁢                                      xe2x80x83                                    ⁢                                      L                    e                    2                                                              ⁢                              xe2x80x83                            +                              xe2x80x83                            ⁢              1                                      ⁢                  
                ⁢                  or          ,                                    (        17        )                                A        ⁢                  xe2x80x83                ≅                  xe2x80x83                ⁢                                            k              e                        ⁢                          xe2x80x83                        ⁢                          B              h              2                                                          k              h                        ⁢                          xe2x80x83                        ⁢                          B              e              2                                                          (        18        )            
Inspection of Eq. (16) leads to the following interesting observations:
The relative spectral centroid shift can be used as a direct strain estimator. We can also observe a direct analogy between the classic definition of strain (Eq. (4)) and the estimator of Eq. (16), which establishes this method as a simple and straight-forward way of estimating the strain.
When the strain is positive (or, compressive), a frequency upshift occurs, i.e.,fc2xe2x88x92fc1 greater than 0. Conversely, a tensile (or, negative) strain results in a frequency downshift, i.e.,fc2xe2x88x92fc1 less than 0. Therefore, the estimator of Eq. (16) provides directly not only the magnitude of the strain but also its sign.
Since constant A is independent of the strain, it will introduce a uniform bias on the resulting elastogram. This should not affect the resulting elastogram, since the latter depicts relative values of strain. The reader should note that the effect of local bandwidth variations is ignored.
The scattering spectrum must be a bandpass and bandlimited spectrum in order to estimate the strain using the centroid estimator. Otherwise, if Be, is infinite and/or if ke is zero, constant A in Eq. (16) will always be zero regardless the strain.
For relatively larger strains, using the more general form of Eq. (3), Eq. (16) becomes                     f        c2            -              f        c1                    f      c1        ≅      A    ⁢          s              1        -        s            
or, by solving for the strain,                     s        ≅                                            f              c2                        -                          f              c1                                                                          (                                  A                  -                  1                                )                            ⁢                              f                c1                                      +                          f              c2                                                          (        19        )            
which is a less straight-forward, but still a direct way of estimating higher strains.
Eq. (16) is reminiscent of the well-known Doppler effect, according to which the ratio of the centroid shift to the center frequency may provide a reliable measure of velocity. According to Eq. (16), in the case of strain, a similar effect also occurs. Also, among others have shown how in broadband Doppler the bandwidth of the resulting spectrum also changes with velocity and that the output RF Doppler spectrum is a frequency-shifted and compressed (or stretched) replica of the transmitted one. Similarly, in the case of strain measurement, that the following expression provides a direct estimation of strain:                                                         B              2                        -                          B              1                                            B            1                          ≅        ks                            (        20        )            
where and are the post- and precompression bandwidths, respectively and k is a constant. So, strain, like velocity, introduces these two coupled effects of centroid shift and bandwidth variation in the power spectrum. Spectral broadening (i.e., in the case of compression) or contraction (i.e., in the case of tension) can introduce a bias in the measurement. A general expression is derived, linking the centroids fc1 and fc2 (before and after compression, respectively), the pure frequency shift (i.e., in case the pre and postcompression spectra are identical, only centered at different frequencies separated by a shift), xcex94f, and a bias term xcex2 denoting the spectral broadening (or, compression) due to strain s:
fc2xe2x88x92fc1=xcex94f+xcex2(s)xe2x80x83xe2x80x83(21) 
where
xcex2(s)=B2(s)xe2x88x92B1 xe2x80x83xe2x80x83(22) 
In order to estimate the strain without the bias associated with spectral broadening, the following equation can be used that results from Eqs. (16), (20), (21) and (22):
xcex94f≅Afc1sxe2x88x92kB1sxe2x80x83xe2x80x83(23)
and solving for strain, the unbiased estimator is given by                     s        ≅                              Δ            ⁢                          xe2x80x83                        ⁢            f                                              Af                              c                1                                      -                          kB              1                                                          (        24        )            
However, Eq. (24) requires a bandwidth estimation and since the bandwidth estimator is not part of this study, we use Eq. (16) as the strain estimator and show a bias with simulations, which is partly due to the previously described bias due to spectral broadening.